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1 his paper covers the design of a dynamic positioning system for the U.S. Coast Ciiiard 
\\’LB "IRIS" class buoy tender. I'he control system design is based upon optimal con- 
trol theory with estimates of position and heading provided by a steady state Kalman 
filter. Sea current estimates are provided by a Kalman filter predictor based upon the 
innovations process. The vessel and dynamic positioning system arc simulated using the 
Dynamic Simulation Language (DSL). 



THESIS DISCLAIMER 



The reader is cautioned that computer programs developed in this research may not 
have been exercised for all cases of interest. While every effort has been made, within 
the time available, to ensure that the programs are free of computational and logic er- 
rors. they cannot be considered validated. Any application of these programs without 
additional verification is at the risk of the user. 
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I. INTRODUCTION TO DYNAMIC POSITIONING SYSTEMS 



A. DEFINITION 

A dynamic positioning system (DPS) is an automatic control system that will drive 
a vessel to a selected position and maintain that position and selected heading using 
computer control of propulsive thrusters. A DPS includes one or more position and 
heading measurement systems and a computer-directed control system. The system 
should keep the vessel within its specified position and heading constraints with optimal 
use of propulsors, reducing the wear on that machinery. 

B. ADVANTAGES AND DISADVANTAGES 

The use of DPS is not new. The first dynamic positioning systems appeared in the 
late 1960 s on coring, cable laying and surface-support ships for underwater work. 1 he 
success of these initial DP vessels spread the application of the systems into ofi'shore oil 
industry, dredging, precision dumping, pipe laying, and (loating hotels. 

1 here are several advantages of dynamic positioning over anchor deployment. It is 
usable in all water depths and not constrained by the length of its anchor chain. In 
contrast to anchoring, there is virtually no set-up time required to begin stationkceping. 
In order to reach and attain a precision position by anchoring, time and personnel in- 
tensive anchoring details must be set. A DPS allows a vessel to attain virtually any po- 
sition and heading automatically. Due to its dynamic nature, a DPS vessel can work in 
close quarters with fixed platforms or anchored vessels. A dynamic positioning system 
is designed to operate in adverse weather conditions 

The disadvantages of DPS include high initial capital investment, high fuel costs 
(over anchoring), high maintenance due to the continuous activation of thrusters, and 
increased manpower requirements for maintenance. 

Dynamic positioning is an elfective means of stationkeeping when water depth is too 
great for anchoring, on-station time will be brief, the vessel is required to be on station 
in all weather conditions, and accuracy in positioning is more important than cost. 

C. EQUIPMENT FOR DYNAMIC POSITIONING SYSTEMS 

The equipment needed for a dynamic positioning system includes a wind sensor 
(velocity and direction), gyrocompass, and a position reference sensor. 
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There are several position sepsors that are employed in DPS designs. They include 
sensors that determine position, velocity and/or acceleration in some conunon reference 
system. The basic position sensor types are taut vi’ire, optical, acoustic beacon, and ra- 
dio systems. , 

The taut wire scheme works by lowering a weight to the sea floor by wire rope. I he 
wire is held in constant tension, usually by a constant-tension winch, and any angle de- 
viations from vertical are corrected by thrusters. 

Optical positioning systems are of either active or passive design and make use of 
laser technology. Triangulation of sextant angles provides position information. How- 
ever, sextants are manually trained by personnel. 

Acoustic beacon is the most common positioning system in the offshore industry. 
This system requires a sea floor beacon and hull-mounted subsurface hydrophones. 

In general, standard radionavigation systems such as Loran-C, OMHGA, and Tansit 
SATNAV do not provide suflicient accuracy for DP systems. 1 he approximate accura- 
cies of these systems are; 

• Loran-C - 460 meters, 

• OMEGA - 500 meters, 

• SATNAV - 25 meters. 

[Ref Ir: p. 8-18] 

Accuracies of those magnitudes are completely unacceptable for DP systems. The 
NAVS'fAR Differential GPS, however, is able to provide accurate position information 
to within 3 meters any^vhere on the globe. In addition, this system allows vessels to 
measure velocity using Doppler to accuracies of about 0.1 meters per second. [Ref 1 : 

p. 8-18). 

For a ship to be dynamically positioned it must use thrusters, sensors, computer, 
and a control and display system. The thrusters must be able to provide fore and aft, 
athwartships, and moment control. Sensors provide physical information including po- 
sition, heading, and wind conditions. The computer uses sensor information and calcu- 
lates required thruster commands that will position the vessel at the reference position 
and heading. Controls and displays are required to allow the operator to monitor and 
control the DP system. 

The purpose of this thesis is to design and simulate a dynamic positioning system 
for a U.S. Coast Guard buoy tender. The dynamic positioning system is designed to be 
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I 'ulK au:ofi'.ated. rerforiiung its functions in van mg sea states, wind, and current con- 
vi'.^ions with Its only operator input being the referenced position. 
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II. BACKGROUND 



The U.S. Coast Guard has begun the acquisition process for seagoing and coastal 
buoy tenders to replace the aging 157 foot WLiM class and the ISO foot WLH class 
tenders. See Figure 1 for a comparison of the old WLB fleet with its successor [Ref 2: 
p.7J. One of the requirements for the replacement vessels is a dynamic positioning sys- 
tem. The dynamically positioned vessel must be able to approach, maneuver, and 
maintain position, within a circle of a 10 meter radius over a fixed point on the earth, 
alongside floating aids to navigation in restricted channels of shallow, 18 foot deep bays, 
estuaries, and rivers. The system is required to maintain position in a 30 knot wind of 
any aspect and 5 knot current within 10 degrees of ship's head [Ref 3 : p. 15 ]. 1 his 

DPS design endeavors to take a modified existing U.S. Coast Guard (USCG) buoy 
tender using a minimum of additional control and measurement inputs in order to reduce 
the potential implementation costs. This initial design is a 180 foot WLB "IRIS" class 
buoy tender. The only modifications to existing vessels are the installation of a stern 
thruster and Global Positioning System (GPS). The GPS provides real-time positioning 
information accurate to within three meters an>"where on the globe [Ref 1 : p. 8-lS ]. 
■fhc GPS is a requirement of the replacement buoy tender fleet. The available controls 
for this design are 200 SUP bow and stern thrusters and 1170 SUP single screw propul- 
sion. System sensor measurements are GPS, gyrocompass, and wind inputs. The mul- 
tiple mission requirements of the buoy tender, as with all Coast Guard cutters, placed 
some constraints on the DPS design from the outset. The buoy tender is occasionally 
called upon to perform some harbor icebreaking duties. It is considered an ice-caixible 
vessel, able to break up to one foot of ice at three knots [Ref. 3 : p. 15 [. Ice operations 
preclude the use of subsurface position sensing devices such as hydrophones and taut 
wire schemes. In order to clarify the reasoning behind some of the system constraints, 
a brief qualitative description of buoy tending operations follows. 

A. BUOY TENDING OPERATIONS 

At the present time, buoy tending operations are at best an inexact science, as is the 
case in many shiphandling routines. The manual ship positioning techniques of experi- 
enced deck officers have been raised to an art form. Put simply, by the use of sextant 
angles or other "local" points of reference, the vessel is slowly aimed at its intended po- 
sition. As the distance to the desired location decreases, the engines are slowed or 
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Figure 1, LB Comparison. (From Ref. 2) 



stopped altogether. As the vessel crosses this desired location, the buoy is released in 
position and the vessel continues on. If, however, the vessel does not cross the desired 
buoy position, another pass must be made. As weather conditions deteriorate, this 
procedure can rapidly increase on station time. A dynamically positioned vessel can be 
prograiruned for a given location and heading. The DP vessel adapts to environmental 
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ciumgcs and nuikes opLiinal use of thrusters to maintain position so it dearly has the 
advantage over the hit-or-niiss scenario just described. 
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III. EQUAl IONS OF MOTION 



1 he niathemniical modeling of the "plant", the waterborne vessel and its environ- 
ment. is divided into several pans: low and high frequency vessel models, wind model, 
and current model. 

A. COORDINATE SYSTEM 

For modelling purposes, two dilfcrent coordinate systems arc used. Earth-lixed co- 
ordinate axes are used to integrate the dynamic equations of motion. A vessel parallel 
set of axes, with the same origin as the earth-fixed axes, is used for computing the forces 
acting on the vessel. See Figure 2 for the relationship of these coordinate axes [Ref. 4; 
p.S53]. The x and y axes represent the vessel in surge and sway and the angle i// is the 
vessel yaw in degrees true. Hence, heading angle is the same for both sets of axes. 
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B. LOW FREQUENCY MOD^L 

The equations of motion for a floating vessel in a seaway arc highly coupled and 
non-linear. Figure 3 is a block diaqram of the low frequency (LF) motions of the vessel 
in surge, sway, and yaw. The low frequency model represents the dynamic response of 
} the ship in "calm" water with wind, current, and propulsive force disturbances. The ef- 
fect of wave action is not considered here as that is modeled by the high frequency model 
below. The coupled non-linear dilTerential equations for surge, sway, and yaw may be 
expressed as follows; 



Xl = -dx\XL~ I - Uc) +Kx + '/£jc (0 

h = 1 'Al -dy I A “ I CVT “ ^'c) + ^wy +^rl + V Ly (2) 

-^’ c ) -d^l^tlh ^ Uk (^) 



where 

• Xl, - LF surge position and velocity 

• Tli A - LF sway position and velocity 

• lA, - LF yaw heading angle and angle rate 

• u^, V, - Current velocities in surge and sway 

• At* Aj - Wind force in surge and sway 

• Fp, Fp - Prop and thruster forces 

• Mp. M, - Wind, thruster, and current moments in yaw 

• Vty '/l'A - Zero mean Gaussian White Noise Processes 

• A. A' A ■ CoelFicients of Drag 

• 7ji^ - Effective Mass Terms 



[Ref 4 : p. 853] 




C. ^ESSEL CHARACTERISTICS 

As prcviou«;ly slated, the modelled vessel is a LSCG WLI3 TRIS " class. T he vessel 
IS considered to ha\e the rollowing characteristics: 

• Length (Overall) - ISO feel 

• Length (betueen perpendiculars) - 1~0 feet 

• Beam-3"y'ce/ 
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• Draft - \2feet 

• Displacement - 943.5 L.Tons 

• m, = 71920 slugs 

• nty = 147615 slugs ’ 

• = l.CilxW slugfi^ 

• = 169 sluglfi 

• = 869 sluglfi 

• d^ = \.\lxW slugfljrad^ 

[Ref. 3 : p. 17] 

The efTective mass (m,, m^„ m^) terms are generated by a computer simulation known 
as the Standard Ship Motion Program (SMP). This program was developed at the 
David W. Taylor Naval Ship Research and Development Center. I'he SMP provides 
predictions of vessel motion in translation (surge and sway) and heading (yaw) in irreg- 
ular seas based upon strip theory. The effective mass is the sum of the vessel mass (the 
mass of water displaced by the vessel) and an "added mass" term. The added mass term 
is predicted by the S.VIP versus varying frequencies of wave encounter, o. Appendix A 
is the printout of the SMP for the 180 foot WLB buoy tender. The added mass and drag 
terms ^lown there are non-dimensional. These terms are a function of the frequency 
of encounter of wave motion, o, which is also non-dimensional. For purposes of this 
thesis, constant dimensionalized values of drag and mass are used, based upon an "av- 
erage" wave frequency of 1.1 radians per second which corresponds to Sea State 5 and 
an average wave height of 5.0 ft. However, the drag terms of the S.MP cannot be di- 
rectly related to the LF equations of motion used here. The drag coeflicients are derived 
directly from Bernoulli's equations for fluid flow on a moving object; 

F=ll2pV^ (4) 

where p is the fluid density and U is the velocity of the fluid relative to the object. 1 he 
drag terms are the product of the Bernoulli force terms with the effective underwater 
area of the vessel. By using the length between perpendiculars, beam, draft, and a block 
coefficient of 0.437, the drag terms are as listed above. Simulations using the lODE 
software package confirm these values. See Appendix B for the simulation results. 'Hie 
vessel is simulated at zero sea state with maximum thrust. Using the effective mass 
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tcr-ns *>• SMP. tr.e drag cccn'r-acv.'s are tliei. 2r:ad un':I li.: k-.ir.vn insvniui: . ’.c/ '.r.c? 
are cr»:aincd. 

D. HIGH FREQUENCY MODEL 

Tiic h:s!. ita^ueucy ves^^c- me Jet rerre?er\:? the hrgh fre^ueiie} 'va^e artiati o: d.e 
vc'^ci Three harnicnic csciUatcrs are modeled m surge, s^vay. and yaw ivuii var.able 
Frc ,uciK>. as (oUc' s 

-'^K ^ ~^H.x ' - ^ 



<h) 



'VI; ere 

• .r,.jr, - HI s. rge pc . .. i:.J • •:!•? ,r> 

• ' - - - HF s '. ay pcs::: :. a:.d ve.c^jy 

• . < . - HE ang'.e ar.d a:.g:e ra*e 

• ■ ‘ - '.y. ■ . - HI' arrga-ar re'.i-cf'ocs :r: surg.. s’^c: . a::i ye v 

• 1. ■ r,, ' . V.. - Zerc ncecL Gi'iss an V. u:e Noise Frorcsses 

iRcl - pSf-i 
E. WIND MODEL 

The w’lid !' nccdel-ed as *.'.;e sum cf siewW vaiying. low mean I're^uc-r.cy ana r^p-di' 
S3r->.?rg i::c ' lre.,uenc} gusts, bod: in v.mJ speed and dj-c-non. Tiie erjuacjons icr the 
: id rncdc.' arc 



'111 

i:i 
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^ds = Vwds 



(13) 



’'V ~ ''V + Wc// ( • 

j where the model terms are defined as follows: 

• w„ - Slowly varying wind velocity 

• - Rapidly varying wind velocity 

• - Slowly varying wind direction 

• w^f - Rapidly varying wind direction 

• ^Kdf - Zero mean Gaussian White Noise Processes 
The wind force terms of equations 1, 2, and 3 are then defined by: 



Kx = + ^v) 


(15) 




(16) 




(17) 


= ^^'ds+'^^df-{'pL + ^h) 


(18) 



where 

• - Wind forces and moment 

• /? - Wind angle of attack 

• F,, Fj, F 3 - Effective wind surface areas 

[Ref 4 : p. 854] 

The F, terms are again derived from Bernoulli's equations for fluid flow. In this case 
the fluid is air and the area is the effective area of wind resistance of the vessel. 

F. CURRENT MODEL 

The current model is described by the following set of differential equations: 



^CE - ^xCE 


(19) 


y'CE = ^yCE 


( 20 ) 
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'I'CE — ^^xj/CE 



( 21 ) 



cos \jy + sin \p 

IV = -XcE sin <A +;^c£ cos (A 
with the terms defined as follows: 

• -^cE-ycE * Current velocities in earth coordinates (North and East respectively) 

• u„ IV - Current velocities in surge and sway 

[Ref 4 ; p. S54] 

1 hese wind, current, high frequency, and low frequency models constitute the vessel 
"plant". Each model represents a different aspect of the plant dynamics. It is clear, 
however, that the only model Avhich may be controlled by the DP system is the low fre- 
quency model. 



( 22 ) 

(23) 
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IV. CONTROL SYSTEM DESIGN 



A. PROPULSORS 

The design of the buoy tender control system is based on several constraints. First, 
there are only three propulsive forces available: the propeller for thrust for\vard and 
astern {±x direction, surge) and the bow and stern thrusters for moment control 
(+ i/^direction, yaw) and athwartships motion (+ y direction, sway). 

The bow and stern thrusters are 200 HP fixed tunnel units. A good rule of thumb 
for the thrust developed by fixed azimuth tunnel thrusters is 30 ftlbfjSllP. These will 
yield approximately 6000 pounds force in the ± y direction. The bow thruster is located 
78.88 feet forward of the center of buoyancy and the stern thruster is 81.12 feet astern 
of it. See Figure 4 for the thruster layout. 

The force and moment equations for the bow and stern thrusters are thus: 



FT=BT+ST 


(24) 


3/1= 78.8857 - 81. 125F 


(25) 



where^ 

• FT - Athwartships force due to thrusters 

• MT - Yaw moment due to thrusters 

• BT - Force of bow thruster 

• ST - Force of stern thruster 

It must be noted that athwartships force, FT, is assumed to be applied at a single point, 
the center of buoyancy. The moment force, MT, is applied about the center of 
buoyancy by the thrusters proportional to their output and the length of the moment 
arm. This length is measured from the center of buoyancy to the thruster. A positive 
moment is considered to be one which will make the ship turn to starboard, a negative 
moment to port. See Figure 4 for the thruster moment dimensions. 

The fixed pitch propeller is driven by a single 1 170 SHP diesel electric power plant. 
It can normally produce upwards of 20,000 pounds force of thrust. The equation for the 
maximum thrust generated by the propeller is: 
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Figuie 4. 1 111 ustcr Forces and Moments 



(«') 

iNssmne a propeller efliciency, = .75 and a maxiinum velocity, = 14.3 kiwis . this 
yields a maximum thrust FP^„= 20,250 fi Ibf. However, due to the natuie of 
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posttionkceping, tlie propeller controller was limited to 12,000 pounds force in the DP 
mode. 

The controller design is based only on the low frequency vessel model. Responding 
to the high frequency wave action of the vessel would cause the controllers to modulate 
in even moderate seas. 

1 he low frequency system was placed in state variable format, 



X = Ax -f Du + Gw 
y = Cx + Du + V 

Ihe non-linear, coupled A and 13 matrices arc: 



1 

-d^\x- u,\liii^ 
0 

0 



0 

1 

0 



- {'fly - 0 -ci^ I .// 1 Ini^ 



1 he state observations matrix, C is: 



The low frequency plant is linearized about 



0 


0 


0 


1///I, 


0 


0 


0 


0 


0 


0 


llniy 


llniy 


0 


0 


0 


0 


78.88 


-81.12 


'1 0 


0 0 


0 o' 


0 0 


1 0 


0 0 


0 0 


0 0 


1 0 



(27) 

(2S) 



(29) 



(30) 



(31) 



( 32 ) 



and the vessel parameters substituted to yield the linear A, B, and G matrices. 
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0 


1 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 



B = 



0 


0 


0 


1.39£-05 


0 


0 


0 


0 


0 


0 


6.11 E- 06 


6.11 E- 06 


0 


0 


0 


0 


78.88 


-81.12 



The plant noise matrix is: 

0 

.0005 

0 

G = 

.0005 

0 

.0005 



The state variable matrix is: 



( 33 ) 



(34) 



(35) 



X 



y 

y 



4 ' 



(36) 
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The control matrix is: 



rp 



I 



u = 



BT 

ST 



(37) 



The Control Toolbox of the software package PC-MATLAB (for MS-DOS Per- 
sonal Computers, Version 3.2-PC June 8, 1987) is then used to obtain the feedback 
gains. Assuming full state feedback is available, the Iqr (linear quadratic regulator) 
function calculates the optimal feedback gain matrix such that the feedback law: 



u = -K^ptX 



(38) 



minimizes the cost function: 



■f'' 



J = I [x^Qx + Ru)di 



where the state weighting matrix, Q is: 



3 £04 
0 
0 
0 
0 
0 



0 

3 £04 
0 
0 
0 
0 



0 

0 

6£04 

0 

0 

0 



0 

0 

0 

6£05 

0 

0 



0 

0 

0 

0 

I £09 
0 



0 

0 

0 

0 

0 

1£09 



and the control weighting matrix, R is: 



R = 



I 0 o' 
0 1 0 
0 0 1 



(39) 



(40) 



(41) 
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The Iqr function returns the following feedback gains: 



} 



173 


4977 


0 


0 


0 


0 


0 


0 


150 


4434 


24962 


249962 


'0 


0 


193 


5701 


-19415 


-19415 



This places the eigenvalues of the closed loop linearized system at; 



-0.035 ±y0.035 
-0.034 ±y0.034 
- 1.0 
-3.6£06 



(42) 



(43) 



Since all the states are not available for feedback, an observer is utilized to provide 
an estimate of the states. A Kalman filter is designed with the aid of the PC-MA'I LAB 
function Iqe (linear quadratic estimator). Recall the system; 

X = Ax Bu Gw (44) 



and 



y = Cx + Du y V (45) 

where the expected values of the plant noise, w, and the measurement noise, v, are 
Zr[u] = ZT[v] = 0 and the covariance matrices for plant and measurement noise are, re- 
spectively. £[wir'] and The covariance matrices are: 



£[ww'] =0.0005 



(46) 



£[n'] 



100 0 0 
0 100 0 
0 0 .0003 



(47) 



The Iqe function returns the Kalman filter gains, K, for the equation: 



X = Ax y Bu y K{y — Cx) (48) 

where x is an estimate of jr. The Kalman gain matrix, K, is: 
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m 



0.0669 0 0 

0.0022 0 0 

0 0.0669 0 

0 0.0022 0 

0 0 1.6069 

0 0 1.2910 



The output of the Kalman filter is the sum of the nonlinear low frequency system, in- 
cluding sampled values of position, heading, and wind velocity and direction, with the 
Kalman updates of position. The sampled wind velocity and direction data is trans- 
formed into ship's coordinates and forces the vessel Kalman filter exactly, within sam- 
pling error, as the true wind forces the low frequency model. The optimal feedback gains 
are then used with the Kalman filter output to create the optimal Kalman filtered 
control: 



The reference error matrix, e, is defined by: 



' x ~ 




V->A 



(50) 



( 51 ) 



The r values here represent the reference positions, heading, and associated velocities. 
The reference positions and heading are the desired location of the vessel. Since the 
desired location is to be maintained, clearly the reference velocities must be zero. 

These optimal controls are then passed to a non-linear filter. This non-linearity is 
a limiter which prevents the thrusters and main propulsion from exceeding the maximum 
DP mode limits (i.e., 12,000 /i>/Cor main propulsion and 6,000 /bf Cor thrusters ). 

It is readily apparent that the control system as designed does not correct for the 
steady state error caused by the introduction of current. A modification to the control 
design to compensate for current offset is imperative. A Kalman filter predictor of 
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current is designed based upon the residual error in the vessel position estimator. That 
is, the diflcrcnce in estimated and measured position and heading is assumed to be due 
to current. This is the "measured" value of the Kalman current filter. Recall the Kalman 
filter equation: 



X = Ax + Bu + K{y — Cx) 



For the current predictor filter: 




( 52 ) 



( 53 ) 



3Y-T 
07 - 0 



0 1 0 0 0 0 
0 0 0 0 0 0 
0 0 0 1 0 0 
0 0 0 0 0 0 
0 0 0 0 0 1 
0 0 0 0 0 0 



c = 



10 0 0 0 0 
0 0 1 0 0 0 



0 0 0 0 1 0 



( 54 ) 



( 55 ) 



( 56 ) 
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The current estimator covariance jnatrices are: 



£[ww'], = 0.00002 



( 57 ) 



5 0 0 



£[vv'],= 0 5 0 

0 0 .005 



(58) 



The current estimation filter gain matrix, K„ is: 



0.0140 0 0 

0.0044 0 0 



0 0.0140 0 

0 0.0044 0 

0 0 6.0000 

0 0 18.200 



(59) 



I he terms and \j/j- are the observations of surge, sway, and yaw including high 

frequency wave action and measurement noise. The estimates of current are converted 
to force terms by Bernoulli's equation and subtracted from the appropriate thruster 
control commands. Figure 5 shows the final system control design in block diagram 
form. ‘ 



22 



Figure 5. Control System De«;ign 



V. COMPUTER SIMULATION 



The equations of motion and control system discussed in previous chapters tvere 
implemented using the Dynamic Simulation Language/VS (DSL). It is a high, level 
continuous simulation language developed by IBM which incorporates VS FORI IU\N 
as a subset. The program is included in Appendix C here as an immediate reference. 
Comments are included as additional documentation. DSL was used here since it is a 
high level language w^hich models the system in continuous time. Although the system 
is modeled as an analog system, DSL allows for the sampling of data as well. This 
sampling was used w'ith the position, heading, and wind measurements, which were up- 
dated at 1 second intervals. Nearly all of the various models require Gaussian noise in- 
puts which are represented by the NORMAL function. This function has inputs of seed 
number, mean, and standard deviation. One final w’ord on DSL is that it has its own 
plotting routine. System outputs presented here were printed using this plotting 
package. 
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VI. RESULTS 



The simulation results oF the dynamically positioned buoy tender are shown here in 
a number of different environments. As discussed previously, the specifications for the 
desired dynamic positioning system are ambitious. To review, the vessel is to attain and 
maintain position to within a 10 meter radius in 30 knots of wind in any aspect and 5 
knots of current within 10 degrees of the bow. The initial test of the simulation is with 
no current, no wind, and at sea state 0. These conditions introduce the least amount 
of disturbance into the system. As can be seen from Figure 6, the vessel holds position 
in surge to within 5 feet. Similarly, Figure 7 and Figure 8 show the vessel holding sway 
position to within 5 feet and a decaying exponential response in yaw angle, with its 
maximum overshoot of 0.05 radians. In each graphical presentation of vessel position 
and heading, the true-measured position and heading are plotted along with the respec- 
tive vessel Kalman filter estimates. A position plot is shown in Figure 9. The position 
plot is presented as an "overhead" look at the vessel track in vessel coordinates, in order 
to obtain this vessel response, a conservative use of propulsive force is used, as witnessed 
in Figure 10 and Figure 11. 

The system is again simulated at sea state 0 with no wind and a current of 3 feet per 
second.' In Figure 12 the system performance is observed in the surge direction. Since 
the vessel estimator is not forced by the current estimator, there is a bias introduced by 
the vessel estimator. This bias could be corrected by modifying the state variable matrix 
with the current velocities in surge and sway such that: 



X 



X — Uc 



y 

y - ^'c 






(60) 



This bias is exactly balanced in the propulsion scheme by the current estimated force. 
Similarly, the system performance in sway. Figure 13, and yaw. Figure 14, is within the 
system constraints in the presence of this current. The position plot. Figure 15, closely 
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resembles its counterpart for the zero current case. However, tlie overshoot in both 
surge ami sway is even more apparent in this plot. The transient performance is not of 
particular importance as the steady state performance is acceptable. Due to the con- 
tinuing nature of the current disturbance, non-zero propulsive forces result in the steadi- 
state as seen in I'igure 16 for main propulsion force and in Figure 17 for bow and stern 
thruster forces. 

The vessel when simulated at a curreiit velocity of 4 feet per second docs not main- 
tain its desired position due to the steady state saturation of the bow and stern thrusters. 
The maximum current for whicli this .system is effective is considered to be 3 feet per 
sccoml. 

Since fully developed seas arc correlated with wind velocity, the final simulation is 
at sea state 5 with current of 3 feet per second and 30 knot winds. With the addition 
of the Kalman predictor for current and the fact that wind is sampled each second and 
is incorporated into the vessel position Kalman filter, no great degradation in system 
performance is anticipated. The surge, sway, and yaw performance of the vessel can be 
seen in Figure 18, Figure 19, and 1 igurc 20, respectively. I he vessel response in surge 
and sway arc well within the specifications. Although thers is no requirement to main- 
tain vessel heading, the steady state yaw response of the vessel is held to about 0.2 radian 
of the reference. The load on the propulsors under these conditions is hca \7 as can be 
seen in Figure 21 and Figure 22. The position plot of the vessel under these conditions 
is seen in Figure 23. The steady state position is still within the vessel specifications. 

Although there is a degradation in system performance with increased environ- 
mental disturbances, the vessel is able to maintain a desired position to within the con- 
straints imposed at the outset. 
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Figure 6. Response in Surge: Sea state 0, No current, No wind. XT is measured 

surge position and XE is estimated surge position. 
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Figure 7. Response in Sway: Sea state 0, No current, No wind. 'i"I is measured 

sway position and YE is estimated sway position. 
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Figure 8. Response in Van: Sea state 0, No current, No wind. PI' is measured 

yaw angle and PE is estimated yaw angle. 
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Figure 9. 



Position Plot: Sea state 0, No wind, No current 



30 




Figure 10. Main Propulsion Force: Sea state 0, No wind, No curreiv 
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Figure 1 1. Do>\ and Stern Thruster Forces: Sea state 0, No wind, No current. 
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Figure 12. Surge Response: Sea state 0, No wind, Current = 3 fps. XT is 

measured surge position and Xli is estimated surge position. 
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Figure 13. Snay Response: Sea state 0, No wind, Current = 3 fps. VT is 

measured sway position and YE is estimated sway position 
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PT VS TIME 



Figure 14. \ Response: Sea slate 0, No wind, Current = 3 fps. PT is meas- 

ured yaw angle and PE is estimated yaw angle. 
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Figure 16. Main Propulsion Force: Sea state 0, No wind, Current = 3 fps 



37 




Figure 17. Bo>> and Stern Thruster Forces : Sea stale 0, No wind. Current = 3 

Ips 
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Figure 18. Surge Response: Sea state 5, Wind = 30 knots, Current = 3 fps. 

XT is measured surge position and XE is estimated surge position. 
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Figure 19. Suay Response: Sea state 5, Wind = 30 knots, Current = 3 fps. YT 

is measured sway position and YE is estimated sway position. 
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1 

Figure 20. Yaw Response: Sea state 5, Wind = 30 knots, Current = 3 fps. PT 

is measured yaw angle and PE is estimated yaw angle. 
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Figure 21. Main Propulsion Force: Sea stale 5, Wind = 30 knots, Cun cm = 3 

(ps 
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Figure 22. Bou and Stern Tliruster Forces: Sea state 5, Wind = 3U knots, Cur- 

rent = 3 fps 
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1 

figure 23. Position Plot: Sea state 5, Wind = 30 knots, Current = 3 I'ps 
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MI. CONCLUSIONS 



This computer-based dynamic positioning system is a viable alternative to manual 
position keeping methods in environmental conditions of up to sea state 5, current up 
to 1.3 knots, and wind up to 30 knots, particularly when precision reciuircments are 
deemed to be of greater importance than fuel cost. In order to meet the requirments for 
the new buoy tender fleet, bow and stern thrusters olMOO SUP should be installed. 
I hcse will allow the DP system to maintain vessel position and heading in the environ- 
ments specified. In addition, variable azimuth thrusters will provide enhanced control 
system performance over fi.xed tunnel thrusters. 
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APPENDIX A. STANDARD SHIP MO l ION PROGRAM 



The added mass and drag cocfricients for the ISO foot WLB "IRIS" class buoy 
lender are dimensionalized forms of the values presented here. 



Ilill. Ill if ii 

tssSii'il" i:;"!!!; s5S:= :J;i 

iii'I.’ itiiii.,. ii; i:i 






rigure 24. SSMP ^ essel Particulars 
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Figure 25. SSMP Added Mass and Damping Coefficients 
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APPENDIX B. lODE SIMULATIONS 



'I'hc following are printouts from lODH programs. lODE (Interactive Ordinary 
Din'ercntial Equations) is an interactive software package that runs on the V.M.C.MS 
mainframe. lODE interrogates the user about the system of dilfcrcntial equations and 
provides the solution in tabular and/or graphical form. 

The lirst program of this appendix is a model of the low Ircqucncy vessel in surge. 
The vessel is modeled at full propulsive force, with no environmental disturbances. Us- 
ing the SMP value for vessel clfcctivc mass and the derived value for the drag coellicicnt, 
the vessel is simulated in the open loop configuration and velocity plotted against time 
in Figure 26 on page 50 'fhc lODE program listing is: 

VARIABLES & INITIAL CONDITIONS: 

X = . 0 
XD = . 0 
T = . 0 

CONSTANTS: 

MX = 71920. 00000 
NX = . 0 

DX = 169. 0000000 
FWX = . 0 
UC = . 0 

U = 20250. 00000 

SPECIAL FUNCTIONS; 

CC = XD-UC 



DERIVATIVES: 

D(X /D(T ) = = 

XD 

D(XD /D(T ) = = 

NX+( U+FWX-DX*ABS ( CC ) *CC ) /MX 

OUTPUTS: 

TITLE: SURGE NO FEEDBACK 
NO TABULATION 
PLOT: X XD 

AGAINST: T 
PLOT: XD 

AGAINST: T AT INTERVAL . 2000000000 
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Figure 26. Surge Motion «itli No Feedback 
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As in the first simulation, thp low frequency sway of the vessel is modeled in the 
open loop and velocity plotted against time in Figure 27 on page 51. The lODE pro- 
gram for sway is: 



VARIABLES & INITIAL COt^DITIONS: 
Y = . 0 
YD = . 0 
T = . 0 

CONSTANTS: 

MY = 147615. 0000 
DY = 869. 0000000 
FWY = .0 
VC = .0 

FT = 12000. 00000 
NY = .0 



SPECIAL FUNCTIONS: 

CC = YD-VC 

DERIVATIVES: 

D(Y /D(T ) = = 

YD 

D(YD /D(T ) = = 

NY+ ( ( FT+FWY - DY* AB S ( CC ) *CC ) /MY ) 

OUTPUTS: 

TITLE: SWAY NO FEEDBACK 
NO TABULATION 
PLOT: YD 

■AGAINST: T AT INTERVAL . 2000000000 

END CALCULATION WHEN T . GE. 150. 000 
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Finally, the yaw motion of the low Ircquency model is simulated in the open loop 
mode and heading angle rate plotted against time in Figure 2S. 
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VARIABLES & INITIAL CONDITIONS: 

P = . 0 
PD = . 0 
T = .0 

CONSTANTS: 

MP = 267000000. 0 
DP = 1170000000. 

NP = .0 
HW = .0 
MC = . 0 

HT = 954000. 0000 
DUM = . 0 

DERIVATIVES: 

D(P /D(T ) = = 

PD 

D(PD /D(T ) = = 

NP+(MT+MW+HC-DP*ABS(PD)*PD)/MP 

OUTPUTS: 

TITLE: PS I DOT 
NO TABULATION 
PLOT: PD 

AGAINST: T AT INTERVAL 

END CALCULATION WHEN T . GE. 150.000 



. 2000000000 
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APPENDIX C PROGRAM LISTING 



'I he ^ollo\^•ing is a listing of the DSL program whicli simulates the vessel, its envi- 
ronment, and the dynamic positioning control system. Comments are included 
throughout to augment the explanations provided in the body of the thesis. 



^DEVELOPED BY W. R. CAIRNS, LCDR, USCG FOR THESIS RESEARCH: 

*"DESIGN AND SIMULATION OF A DYNAMIC POSITIONING SYSTEM FOR A U. S. 
*COAST GUARD BUOY TENDER" 

TITLE: DPSIM5 

mils PROGRAM SIMULATES A WLB "IRIS" CLASS USCG BUOY TENDER SUBJECT TO 
*WIND, CURRENT, AND WAVE FORCES AND CONTROLLED BY A DYNAMIC POSITIONING 
*SYSTEM 

*THE PROGRAM PARAMETERS (PARAM) ARE DEFINED AS: 

* F1,F2,F3 ARE THE VESSEL EFFECTIVE WIND AREAS IN SURGE, SWAY, AND YAW 

* AI,A2 ARE THE VESSEL WIND VELOCITY AND DIRECTION FORCE COEFFICIENTS 

* RX,RY,RP ARE THE REFERENCE POSITIONS IN SURGE, SWAY AND YAW 
PARAM Fl=2. 294,F2=I1. 16,F3=474. 3,. . . 

A1=0. 001,A2=0. 001,. . . 

W=00. 8,. . . 

IWVS=30. 0,. . . 

IWDS=0. 0,. . . 

XCI=0. 0,. . . 

YCI=0. 0,. . . 

PCI=0. 0,. . . 

RX=]00. 0,. . . 

RY=100. 0,. . . 

RP=1. 00 



*THE CONSTANTS (CONST) ARE DEFINED: 

* X0,XD0 INITIAL CONDITIONS FOR SURGE VELOCITY AND ACCELERATION 

* Y0,YD0 INITIAL CONDITIONS FOR SWAY VELOCITY AND ACCELERATION 

* P0,PD0 INITIAL CONDITIONS FOR YAW VELOCITY AND ACCELERATION 

* MX,DX MASS AND DRAG COEFFICIENTS IN SURGE 

* MY,DY MASS AND DRAG COEFFICIENTS IN SWAY 

* MP,DP MASS AND DRAG COEFFICIENTS IN YAW 

* NX,UX,SX SEED, MEAN, STD DEVIATION FOR LF PLANT NOISE IN SURGE 

* NY,UY,SY SEED, MEAN, STD DEVIATION FOR LF PLANT NOISE IN SWAY 

* NP,UP,SP SEED, MEAN, STD DEVIATION FOR LF PLANT NOISE IN YAW 

* HNX,HUX,HSX SEED, MEAN, STD DEVIATION FOR HF PLANT NOISE IN SURGE 

* HNY,HUY,HSY SEED , MEAN , STD DEVIATION FOR HF PLANT NOISE IN SWAY 

* HNP,HUP,HSP SEED, MEAN, STD DEVIATION FOR HF PLANT NOISE IN YAW 

* NtYX,UW'X,SWX SEED, MEAN, STD DEVIATION FOR HF FREQ NOISE IN SURGE 

* NWY,UWY,SWY SEED, MEAN, STD DEVIATION FOR HF FREQ NOISE IN SWAY 

* Nt^P,UWP,SWP SEED, MEAN, STD DEVIATION FOR HF FREQ NOISE IN YAW 

* MNX,UMX,SMX SEED, MEAN, STD DEVIATION FOR LF MEAS NOISE IN SURGE 

MNY,UMY,SMY SEED , MEAN , STD DEVIATION FOR LF MEAS NOISE IN SWAY 

* MNP,UMP,SMP SEED, MEAN, STD DEVIATION FOR LF MEAS NOISE IN YAW 

* CNX,CUX,CSX SEED, MEAN, STD DEVIATION FOR CURRENT NOISE IN SURGE 

* CNY,CUY,CSY SEED, MEAN, STD DEVIATION FOR CURRENT NOISE IN SWAY 

* CNP,CUP,CSP SEED, MEAN, STD DEVIATION FOR CURRENT NOISE IN YAW 

* NGX,MCX,SCX SEED, MEAN, STD DEVIATION FOR CURRENT NOISE IN SURGE( KALMAN) 
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* NCY,MCY,SCY SEED, MEAN, STD DEVIATION FOR CURRENT NOISE IN SWAY (KALMAN) 

* NCF,MCP,SCP SEED, MEAN, STD DEVIATION FOR CURRENT NOISE IN YAW (KALMAN) 

* M>VS,UWVS,SWVS SEED, MEAN, SD FOR SLOW WIND VELOCITY NOISE 

* M>VF,UVVF,SWVF SEED, MEAN, SD FOR FAST WIND VELOCITY NOISE 

* NWDS,UWDS,SWDS SEED, MEAN, SD FOR SLOW WIND DIRECTION NOISE 

* NVDF,UWDF,SWDF SEED, MEAN, SD FOR FAST WIND DIRECTION NOISE 

* FPR.FPL PROPULSION LIMITER (+/-) 

* BTR,BTL BOW THRUSTER LIMITER (+/-) 

* STR,STL STERN THRUSTER LIMITER (+/-) 

* EXD,EXU DEAD SPACE IN SURGE 

* EYD.EYU DEAD SPACE IN SWAY 

* EFD,EPU DEAD SPACE IN YAW 

* RXD,RYD,RPD REFERENCE VELOCITIES IN SURGE, SWAY, YAW 

* K VESSEL KALMAN FILTER GAINS SUBSCRIPTED ACCORDING TO MATRIX 

* C CURRENT ESTIMATION FILTER GAINS SUBSCRIPTED ACCORDING TO MATRIX 

* G OPTIMAL FEEDBACK GAINS SUBSCRIPTED ACCORDING TO MATRIX 

CONST X0=0. O.XDO=0. 0 

Y0=0. 0,YD0=0. 0,. . . 

P0=0. 0,PDO=0. 0,. . . 

MX=71920. 0,DX=169. 00000,. . . 

MY=1^7615. 0,DY=869. 0000,. . . 

MP=2. 67D+08,DP=1. 17D+09,. . . 

NX=19. 0,UX=0. 0,SX=0. 0005,. . . 

NY=21. 0,UY=0. 0,SY=0. 0005,. . . 

NF=23. 0,UP=0. 0,SP=0. 00001,. . . 

HNX=19. 0,HUX=0. 0,HSX=3. 0,. . . 

HNY=21. 0,HUY=0. 0,HSY=3. 0,. . . 

HNP=23. 0,HUP=0. 0,HSP=0. 175,. . . 

N'WX=19. 0,UVX=0. 0,SWX=0. 005 , . . . 

M7Y=21. 0,UVY=0. 0,SWY=0. 005,. . . 

N¥P=23. 0,UWT=0. 0,SWP=0. 000175,. . . 

MNX=19. 0.1’1X=0. 0,SMX=3. 0000,. . . 

MNY=21. 0,UMY=0. 0,SMY=3. 0000,. . . 

MNP=23. 0,UMP=0. 0,SMP=0. 0175,. . . 

CNX=.19. 0,CrX=0. 0,CSX=0. 00050,. . . 

CNY=21. 0,CUY=0. 0,CSY=0. 00050,. . . 

CNP=23. 0,CUP=0. 0,CSP=0. 00001,. . . 

NCX=19. 0,MCX=0. 0,SCX=0. 0005,. . . 

NCY=2 1.0, MCY=0. 0 , SCY=0. 0005 , . . . 

NCP=23. 0,MCP=0. 0,SCP=0. 00001,. . . 

N'WVS=19. 0,U\7VS=0. 0,SWVS=0. 0005 

N'WVr=21. 0,UWVF=0. 0,SWVF=0. 0005,. . . 

N'WDS=23. 0,UWDS=0. 0,SWDS=0. 00001,. . . 

NYvDF=2 3. 0,UWDF=0. 0,SWDF=0. 00001,. . . 

V.NN=15. 0,VNU=0. 0,VNS=0. 0005,. . . 

FFR=12000. 0,FPL=-12000. 0,. . . 

BTR=6000. 0 , BTL=-6000. 0 , . . . 

STR=6000. 0,STL=-6000. 0,. . . 

RXD=0. 0,RYD=0. 0,RPD=0. 0,. . . 

K11=0. 0669,K12=0. 0,K13=0. 0,. . . 

K21=0. 0022,K22=0. 0,K23=0. 0,. . . 

K31=0. 0,K32=0. 0669,K33=0.0,. . . 

K41=0. 0,K42=0. 0022,K43=0. 0,. . . 

K51=0. 0,K52=0. 0,K53=1. 6069,. . . 

K61=0. 0,K62=0. 0,K63=1. 2910,. . . 

Cl 1=0. 014,C12=0. 0,C13=0. 0,. . . 
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C21=0. 0044,C22=0. 0,C23=0. 0,. . . 

C31=0. 0,C32=0. 014,C33=0. 0,. . . 

C41=0. 0,C42=0. 0044,C43=0. 0,. . . 

C51=0. 0,C52=0. 0,C53=6. 0000,. . . 

C61=0. 0,062=0.0,063=18. 200,. . . 

011=193. 0,012=5701. 00 ,013=0. 0 ,014=0. 0 ,015=0. 0,016=0. 0,. . . 

021=0. 0,022=0. 0,023=150. 0,024=4434. 0 ,025=24962. 0 ,026=24962. 0 , . . . 
031=0. 0,032=0. 0,033=193. 0,034=5701. 0,035=-19415. 000,036=-19415. 00 

* 011=78. 50,012=9200. 00 ,013=0. 0 ,014=0. 0 , 015=0. 0 ,016=0. 0 ,. . . 

* 021=0. 0,022=0. 0,023=42. 00,024=7200. 0,025=5854. 0.026=8200. 0,. . . 

* 031=0. 0,032=0. 0,033=38. 0,034=6900. 0,035=-4553. 000 ,036=-7500. 00 

* 011=77. 46,012=3335. 50 ,0 13=0. 0 ,014=0. 0 ,015=0. 0 , 016=0. 0 , . . . 

* 021=0. 0,022=0. 0,023=47. 556,024=2483. 80,025=249. 62,026=136. 7,. . . 

* 031=0. 0,032=0. 0,033=61. 14,034=3193. 40,035=-194. 150 ,036=-106. 34 
DYNAMIC 

*L0W FREQUENCY MODEL PLANT NOISE IN SURGE, SWAY, AND YAW 
NOIX=NORMAL(NX,UX,SX) 

NOIY=NORMAL(NY,UY,SY) 

NOIP=NORMAL(NP,UP,SP) 

*HIGH FREQUENCY MODEL PLANT NOISE IN SURGE, SWAY, AND YAW 
HNOIX=NORMAL( HNX , HUX , HSX ) 

HNOIY=NORMAL( HNY , HUY , HSY) 

HNOIP=NORMAL( HNP ,HUP ,HSP) 

’•'MEASUREMENT NOISE IN SURGE, SWAY, AND YAW 
MNOIX=NORMAL( MNX , UMX , SMX) 

MNOIY=NORMAL( MNY ,UMY , SMY) 

MNOIP=NORMAL( MNP , UMP , SMP) 

’■'CURRENT NOISE IN SURGE, SWAY, AND YAW 
CXNOI=NORMAL( CNX , CUX , CSX) 

CYNOI=NORMAL( CNY , CUY , CSY) 

CPNOI=NORMAL( CNP , CUP , CSP) 

’■'KALMAN FILTER CURRENT NOISE IN SURGE, SWAY, AND YAW 
KCX=NORMAL( NCX , MCX , SCX ) 

KCY=NORMAL( NC Y , MCY , SCY ) 

KCP=NORMAL(NCP,MCP,SCP) 

’■'WIND NOISE (SLOW/FAST VARYING) ( VELOCITY AND DIRECTION) 

WVNS=NORMAL( Nl^VS , U\v’ VS , SWVS ) 

WVNF=NORMAL( NWVF , U\^VF , S WVF ) 

WDNS=NORMAL( NWDS , UVv’DS , SWDS ) 

WDNF=NORMAL( NWDF , UWDF , SWDF ) 

’^ESTIMATOR ERROR - INNOVATIONS IN SURGE, SWAY, AND YAW 
INX=XT-XE 
INY=YT-YE 
INP=PT-PE 

’■'CURRENT CORRECTION (AS APPLIED TO LF VESSEL MODEL) 

Q1=XD-UC 

Q2=YD-VC 

’''THRUSTER FORCE = BOW THRUSTER + STERN THRUSTER 
FT=BT+ST 

’''BOW THRUSTER CONTROL 

BTI= G21’''EX+G22’’'EXXD+G23’''EY+G24’^'EYYD+G25’’'EP+G26’’'EPPD 
BTII=BTI-(DY’''(VCH)’''ABS(VCH)+FWYS) 

’V BTII=BTI-FWYS/2. 0 

BT=LI M IT( BTL , BTR , BT I I ) 

’■'STERN THRUSTER CONTROL 

STI= G31’''EX+G32’’'EXXD+G33’^'EY+G34’''EYYD+G35’''EP+G36’''EPPD 
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STII=STI-(DY*(VCH)*ABS(VCH)+FWYS) 

* STII=STI-FWYS/2. 0 
ST=LIMIT( STL,STR,STII) 

*PROP CONTROL 

FPI= G11*EX+G12*EXXD+G13*EY+G14*EYYD+G15*EP+G16*EPPD 
FPI I=FPI -( DX*UC*ABS( UCH)+0. 75*FWXS) 

* FPII=FPI-FWXS 
FP=LIMIT(FPL,FPR,FPII) 

^MOMENT CONTROL 

MT=78. 88*BT-81. 12*ST 

SAMPLE 

’•■-TOTAL VESSEL MOTION MODEL - MEASUREMENTS IN SURGE, SWAY, AND YAW 
XT=X+XH+MNOIX 
YT=Y+YH+MNOIY 
PT=P+PH+MNOIP 

*SAMPLED WIND FORCES AND MOMENT IN SURGE, SWAY, AND YAW 
FWXS=FWX 
FWYS=FWY 
M\v'S=MW 
DERIVATIVE 

*LOW FREQUENCY VESSEL MODEL 
*LF SURGE 

X=INTGRL(XO,XD) 

XD=INTGRL(XDO,XDD) 

XDD=NO I X+ ( FWX+FP - DX*Q 1*AB S ( Q 1 ) +M Y*Q 2*PD+MX’-'VC*PD ) / MX 
*LF SWAY 

Y=INTGRL(YO,YD) 

YD=INTGRL(YDO,YDD) 

YDD=NOI Y+( FWY+FT-DY*Q2*ABS ( Q2 )+MY*UC*PD+MX*Q 1*PD ) /MY 
*LF YAW 

F=INTGRL(PO,PD) 

PD=INTGRL(PDO,PDD) 

PDD=NO I P+ C - ( M Y - MX ) *Q 1 *Q2 - DP*AB S C PD ) *PD+MC+MW+MT ) / MP 
*11IGH FREQUENCY VESSEL MODEL 
*HF SURGE 

X11=INTGRL(0. 0,XHD) 

XHD=INTGRL(0. O.XHDD) 

XHDD= - ( WX** 2 ) *XH+HNOI X - . 2*XHD 
*HF SWAY 

YH=INTGRL(0. 0,YHD) 

YHD=INTGRL(0. 0,YHDD) 

YHDD= - ( W Y** 2 ) *YH+HNO I Y - . 2* YHD 
*HF YAW 

PH=INTGRL(0. 0,PHD) 

PHD=INTGRL(0. 0,PHDD) 

PHDD= - ( WP**2 ) *PH+HNO IP-. 2*PUD 
*HF WAVE FREQUENCY VARIATION MODEL 
*FREQUENCY IN SURGE 

WX=INTGRL(W,WXD) 

WXD=NORMAL( NWX , UWX , SWX ) 

^FREQUENCY IN SWAY 

WY=INTGRL(W,WYD) 

WYD=NORMAL( NWY , UWY , SWY ) 

*FREQUENCY IN YAW 

WP=INTGRL(W,WPD) 
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WPD=NORtlAL( WP , UWP , SWP ) 

*CURRENT MODEL IN EARTH COORDINATES 
XCED=INTGRL(XCI ,XCEDD) 

XCEDD=CXNOI 
YCED=INTGRL(YCI ,YCEDD) 

YCEDD=CYNOI 

PCED=INTGRL( PCI ,PCEDD) 

PCEDD=CPNOI 

*CURRENT MODEL IN SURGE, SWAY, AND YAW 
UC=XCED*DCOS(PT)+YCED*DSIN(PT) 

VC=-XCED*DSIN( PT) +YCED*DCOS( PT) 

MC=PCED 
*WIND MODEL 

*SLOWLY VARYING WIND VELOCITY 
WVS=INTGRL( IWVS ,WVSD) 

VVSD=WVNS 

^RAPIDLY VARYING WIND VELOCITY 
WVF=INTGRL(0. 0,WVFD) 

WVFD=A1*WVF+WVNF 
*SLOWLY VARYING WIND DIRECTION 
WDS=INTGRL(IWDS,WDSD) 

WDSD=WDNS 

*RAPIDLY VARYING WIND DIRECTION 
WDF=INTGRL(0. 0,WDFD) 

WDFD=A2*WDF+WDNF 
*WIND ANGLE OF ATTACK 
BETA=WDS+WDF-PT 

*WIND FORCE EQUATIONS IN SURGE, SWAY, AND YAW 
FWX= -F l*DCOS ( BETA ) *( WVS+WVF ) **2 
F W Y= - F 2''' D S I N ( B E TA ) W V S+W VF ) **2 
M\v''=F 3*D S I N ( BETA ) * C W VS+WVF ) ** 2 
*KALMAN FILTER FOR VESSEL 
^VESSEL KALMAN IN SURGE 
EX=RX-XE 
EXXD=RXD-XED 

AXE=(-DX-'^(XED-0. )*ABS(XED-0. )+MY*(YED-0. )*PED+MX*0. *PED)/MX 

BXE=(FP+FWXS)/MX 

XE=INTGRL(0. 0,XED) 

XED=INTGRL(0. 0,XEDD)+K11*(XT-XE)+K12*(YT-YE)+K13*(PT-PE) 
XEDD=AXE+BXE+K2 1* ( XT - XE ) +K2 2*( YT - YE ) +K2 3* ( PT - PE ) 

^VESSEL KALMAN IN SWAY 
EY=RY-YE 
EYYD=RYD-YED 

AYE=(-DY*(YED-0. )*ABS(YED-0. )+MX*(XED-0. )*PED+MY*0. *PED)/MY 

BYE=(FT+FWYS)/MY 

YE=INTGRL(0. 0,YED) 

YED=INTGRL(0. 0,YEDD)+K31*(XT-XE)+K32*(YT-YE)+K33*(PT-PE) 
YEDD=AYE+BYE+K41*(XT-XE)+KA2*(YT-YE)+K43*CPT-PE) 

*VESSEL KALMAN IN YAW 
EP=RP-PE 
EPPD=RPD-PED 

APE=(-(MY-MX)*(XED-0. )*(YED-0. ) -DP*PED*ABS(PED) )/MP 
BPE=(IWS+MT+0. )/MP 
PE=INTGRL(0. 0,PED) 

PED=INTGRL(0. 0,PEDD)+K51*(XT-XE)+K52*(YT-YE)+K53*(PT-PE) 
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PEDD=APE+BPE+K6 1*( XT-XE ) +K62*( YT-YEi+K63*( PT-PE ) 

^FILTER PREDICTION OF CURRENT OFFSET 
^CURRENT PREDICTION IN SURGE 
CX=INX-XCH 
XCH=INTGRL(0. 0,UCH) 

UCH=INTGRL(0. 0,UCHD)+(C11*CX) 

UCHD=(C21*INX) 

^CURRENT PREDICTION IN SWAY 
CY=INY-YCH 
YCH=INTGRL(0. 0,VCH) 

VCH=INTGRL(0. 0,VCHD)+(C32*CY) 

VCHD=(C42*INY) 

*CURRENT PREDICTION IN YAW 
CP=INP-PCH 
PCH=INTGRL(0. 0,PCHD) 

PCHD=INTGRL(0. 0,PCHDD)+(C53*CP) 

PCHDD=(C63*CP) 

* TRANSFORM BACK TO VESSEL COORDINATES 
MCH=PCHD 

CONTROL FINT1M=800. 0,DELT=. 10,DELS=1. 0 
METHOD RKSFX 

SAVE 10. 00,X,EX,FP,Y,EY,BT,ST,MT,P,EP.XT,YT,PT,XE,YE.PE,XH,YH,PH,. . . 

FWX , FWY , M\v . XDD , XHDD , MC , UC , VC , UCH , VCH , I NX , I NY , INP , XED , YED , PED 
*RINT UC, UCH, VC, VCH 

GRAPH (G1,DE=TEK618) TIME(UN=' SEC ' ) ,XT(UN=' FT' ) ,XE 
LABEL (Gl) XT VS TIME 

GRAPH (G2,DE=TEK618) TIME(UN=' SEC' ) ,YT(UN='FT' ) ,YE 
LABEL (G2) YT VS TIME 

GRAPH (G3,DE=TEK618) TIME(UN=’ SEC' ) ,PT(UN='RAD' ) ,PE 
LABEL (G3) PT VS TIME 

GRAPH (G4,DE=TEK618) TIME(UN=' SEC' ) ,FP(UN=' LBF' ) 

LABEL (G4) PROP CONTROL 

GRAPH (G5,DE=TEK618) TIME, BT, ST 

LABEL (G5) THRUSTER CONTROL 

GRAPH (G6,DE=TEK618) Y(UN='FT’ ) ,X(UN='FT' ) 

LABEL (G6) POSITION X VS Y 

GRAPH CG7,DE=TEK618) XE,XED 

LABEL (G7) PHASE PLANE 

*GRAPH (G7,DE=TEK618) TIME, UC, UCH 

*LABEL (G7) POSITION UC.UCH VS TIME 

*GRAPH (G8,DE=TEK618) TIME, VC, VCH 

*LABEL (G8) POSITION VC, VCH VS TIME 

*GRAPH (G9,DE=TEK618) TIME,INX,INY 

*LABEL (G9) ERROR RESIDUALS 

*GRAPH (G8,DE=TEK618) TIME,Y,YH 

*LABEL (G8) HIGH AND LOW FREQUENCY SWAY 

*GRAPH (G9,DE=TEK618) TIME, P, PH 

*LABEL (G9) HIGH AND LOW FREQUENCY YAW 

END 

STOP 
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